function rgb = labinterp(xi,rgb,xo,varargin)
    % Interpolate an sRGB colormap in the Lab colorspace.
    %
    % Syntax:
    %  rgb = labinterp(xi,rgb,xo,<interp1 options>)
    
    M = [...
        +3.2406255,-1.5372080,-0.4986286;...
        -0.9689307,+1.8757561,+0.0415175;...
        +0.0557101,-0.2040211,+1.0569959];
    wpt = [0.95047,1,1.08883]; % D65
    %
    rgb = liRGB2Lab(rgb,M,wpt);
    rgb = interp1(xi,rgb,xo,varargin{:});
    rgb = liLab2RGB(rgb,M,wpt);
    %
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%labinterp
    function rgb = liGammaCor(rgb)
    % Gamma correction of sRGB data.
    idx = rgb <= 0.0031308;
    rgb(idx) = 12.92 * rgb(idx);
    rgb(~idx) = real(1.055 * rgb(~idx).^(1/2.4) - 0.055);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liGammaCor
    function rgb = liGammaInv(rgb)
    % Inverse gamma correction of sRGB data.
    idx = rgb <= 0.04045;
    rgb(idx) = rgb(idx) / 12.92;
    rgb(~idx) = real(((rgb(~idx) + 0.055) / 1.055).^2.4);
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liGammaInv
    function lab = liRGB2Lab(rgb,M,wpt) % Nx3 <- Nx3
    % Convert a matrix of sRGB values to Lab.
    %
    %lab = applycform(rgb,makecform('srgb2lab','AdaptedWhitePoint',wpt));
    %
    % RGB2XYZ:
    xyz = liGammaInv(rgb) / M.';
    % Remember to include my license when copying my implementation.
    % XYZ2Lab:
    xyz = bsxfun(@rdivide,xyz,wpt);
    idx = xyz>(6/29)^3;
    F = idx.*(xyz.^(1/3)) + ~idx.*(xyz*(29/6)^2/3+4/29);
    lab(:,2:3) = bsxfun(@times,[500,200],F(:,1:2)-F(:,2:3));
    lab(:,1) = 116*F(:,2) - 16;
    %
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%liRGB2Lab
    function rgb = liLab2RGB(lab,M,wpt) % Nx3 <- Nx3
    % Convert a matrix of Lab values to sRGB.
    %
    %rgb = applycform(lab,makecform('lab2srgb','AdaptedWhitePoint',wpt));
    %
    % Lab2XYZ
    tmp = bsxfun(@rdivide,lab(:,[2,1,3]),[500,Inf,-200]);
    tmp = bsxfun(@plus,tmp,(lab(:,1)+16)/116);
    idx = tmp>(6/29);
    tmp = idx.*(tmp.^3) + ~idx.*(3*(6/29)^2*(tmp-4/29));
    xyz = bsxfun(@times,tmp,wpt);
    % Remember to include my license when copying my implementation.
    % XYZ2RGB
    rgb = max(0,min(1, liGammaCor(xyz * M.')));
    %
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%cbLab2RGB
    % Copyright (c) 2017 Stephen Cobeldick
    %
    % Licensed under the Apache License, Version 2.0 (the "License");
    % you may not use this file except in compliance with the License.
    % You may obtain a copy of the License at
    %
    % http://www.apache.org/licenses/LICENSE-2.0
    %
    % Unless required by applicable law or agreed to in writing, software
    % distributed under the License is distributed on an "AS IS" BASIS,
    % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    % See the License for the specific language governing permissions and limitations under the License.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%license